#include"mbpt.h"
#include<iomanip>
#include<cmath>
#include<algorithm>

using std::abs;
using std::min;
using std::max;
using std::setw;

using Eigen::RealSchur;

MBPT::DataType MBPT::cross_coupled1(int a,int b,int c,int d,int J) const
{
  DataType temp=0;
  int jja=HFOrbs[a].jj;
  int jjb=HFOrbs[b].jj;
  int jjc=HFOrbs[c].jj;
  int jjd=HFOrbs[d].jj;
  int j_min=max(abs(jja-jjb),abs(jjc-jjd))/2;
  int j_max=min(jja+jjb,jjc+jjd)/2;
  if( (HFOrbs[a].tz+HFOrbs[b].tz)!=(HFOrbs[c].tz+HFOrbs[d].tz) ) return 0;
  if( (HFOrbs[a].l+HFOrbs[b].l + HFOrbs[c].l+HFOrbs[d].l)%2 ) return 0;
  if(J<abs(jja-jjc)/2) return 0;
  if(J<abs(jjb-jjd)/2) return 0;
  if(J>(jja+jjc)/2) return 0;
  if(J>(jjb+jjd)/2) return 0;
  
  for(int j=j_min;j<=j_max;j++)
    {
      temp+=(j*2+1)*phase((jjb+jjc)/2+j+J)*gsl_sf_coupling_6j(jjc,jja,2*J,jjb,jjd,2*j)*get2B(a,b,c,d,j);      
    }
  return temp;
}

MBPT::DataType MBPT::cross_coupled2(int a,int b,int c,int d,int J) const
{
  DataType temp=0;
  int jja=HFOrbs[a].jj;
  int jjb=HFOrbs[b].jj;
  int jjc=HFOrbs[c].jj;
  int jjd=HFOrbs[d].jj;
  int j_min=max(abs(jja-jjb),abs(jjc-jjd))/2;
  int j_max=min(jja+jjb,jjc+jjd)/2;
  if( (HFOrbs[a].tz+HFOrbs[b].tz)!=(HFOrbs[c].tz+HFOrbs[d].tz) ) return 0;
  if( (HFOrbs[a].l+HFOrbs[b].l + HFOrbs[c].l+HFOrbs[d].l)%2 ) return 0;
  if(J<abs(jja-jjd)/2) return 0;
  if(J<abs(jjb-jjc)/2) return 0;
  if(J>(jja+jjd)/2) return 0;
  if(J>(jjb+jjc)/2) return 0;

  for(int j=j_min;j<=j_max;j++)
    {
      temp-=(j*2+1)*phase((jjb+jjc)/2+J)*gsl_sf_coupling_6j(jjc,jjb,2*J,jja,jjd,2*j)*get2B(a,b,c,d,j);      
    }
  return temp;
}


void MBPT::calOccuNum()
{
  OccuNum.resize(totalOrbitals);
  for(int i=0;i<totalOrbitals;i++)
    {
      OccuNum[i]=0;
    }
  DataType norm=1;//norm of two order perturbation wave function.
  DataType temp1=0,temp2=0,temp3=0;
  for(int h1=0;h1<=FermiSurface;h1++)
    {
      for(int h2=h1;h2<=FermiSurface;h2++)
	{

	  int jj_min=abs(HFOrbs[h1].jj-HFOrbs[h2].jj);
	  int jj_max=HFOrbs[h1].jj+HFOrbs[h2].jj;
	  int Tz=HFOrbs[h1].tz+HFOrbs[h2].tz;
	  int par=((HFOrbs[h1].l+HFOrbs[h2].l)%2);
	  DataType e=HFOrbs[h1].e+HFOrbs[h2].e;
	  for(int jj=jj_min;jj<=jj_max;jj+=2)
	    {
	      for(int p1=FermiSurface+1;p1<totalOrbitals;p1++)
		for(int p2=p1;p2<totalOrbitals;p2++)
		  {

		    if((HFOrbs[p1].tz+HFOrbs[p2].tz)!=Tz) continue;
		    if((HFOrbs[p1].l+HFOrbs[p2].l-par)%2) continue;
		    if(abs(HFOrbs[p1].jj-HFOrbs[p2].jj)>jj) continue;
		    if((HFOrbs[p1].jj+HFOrbs[p2].jj)<jj) continue;
		    DataType temp=(jj+1.)*get2B(p1,p2,h1,h2,jj/2)*get2B(p1,p2,h1,h2,jj/2)/pow( (e-HFOrbs[p1].e-HFOrbs[p2].e),2);
		    if(p1==p2&&h1==h2)
		      {
			temp3+=temp;
		      }
		    else if(p1==p2||h1==h2)
		      {
			temp2+=temp;
		      }
		    else
		      {
			temp1+=temp;
		      }
		  }
	    }
	}
    }
  temp2*=0.5;
  temp3*=0.25;
  temp1+=temp2;
  temp1+=temp3;
  norm+=temp1;

  for(int h1=0;h1<=FermiSurface;h1++)
    {
      for(int h2=0;h2<=FermiSurface;h2++)
	{
	  int jj_min=abs(HFOrbs[h1].jj-HFOrbs[h2].jj);
	  int jj_max=HFOrbs[h1].jj+HFOrbs[h2].jj;
	  int Tz=HFOrbs[h1].tz+HFOrbs[h2].tz;
	  int par=((HFOrbs[h1].l+HFOrbs[h2].l)%2);
	  DataType e=HFOrbs[h1].e+HFOrbs[h2].e;
	  for(int jj=jj_min;jj<=jj_max;jj+=2)
	    {
	      for(int p1=FermiSurface+1;p1<totalOrbitals;p1++)
		for(int p2=p1;p2<totalOrbitals;p2++)
		  {

		    if((HFOrbs[p1].tz+HFOrbs[p2].tz)!=Tz) continue;
		    if((HFOrbs[p1].l+HFOrbs[p2].l-par)%2) continue;
		    if(abs(HFOrbs[p1].jj-HFOrbs[p2].jj)>jj) continue;
		    if((HFOrbs[p1].jj+HFOrbs[p2].jj)<jj) continue;
		    DataType temp=-(jj+1.)*get2B(p1,p2,h1,h2,jj/2)*get2B(p1,p2,h1,h2,jj/2)/pow( (e-HFOrbs[p1].e-HFOrbs[p2].e),2);
		    if(p1==p2)
		      {
			OccuNum[h1]+=0.5*temp;
		      }
		    else
		      {
			OccuNum[h1]+=temp;
		      }
		  }
	    }
	}
    }


  for(int h1=0;h1<=FermiSurface;h1++)
    {
      for(int h2=h1;h2<=FermiSurface;h2++)
	{
	  int jj_min=abs(HFOrbs[h1].jj-HFOrbs[h2].jj);
	  int jj_max=HFOrbs[h1].jj+HFOrbs[h2].jj;
	  int Tz=HFOrbs[h1].tz+HFOrbs[h2].tz;
	  int par=((HFOrbs[h1].l+HFOrbs[h2].l)%2);
	  DataType e=HFOrbs[h1].e+HFOrbs[h2].e;
	  for(int jj=jj_min;jj<=jj_max;jj+=2)
	    {
	      for(int p1=FermiSurface+1;p1<totalOrbitals;p1++)
		for(int p2=FermiSurface+1;p2<totalOrbitals;p2++)
		  {

		    if((HFOrbs[p1].tz+HFOrbs[p2].tz)!=Tz) continue;
		    if((HFOrbs[p1].l+HFOrbs[p2].l-par)%2) continue;
		    if(abs(HFOrbs[p1].jj-HFOrbs[p2].jj)>jj) continue;
		    if((HFOrbs[p1].jj+HFOrbs[p2].jj)<jj) continue;
		    DataType temp=(jj+1.)*get2B(p1,p2,h1,h2,jj/2)*get2B(p1,p2,h1,h2,jj/2)/pow( (e-HFOrbs[p1].e-HFOrbs[p2].e),2);
		    if(h1==h2)
		      {
			OccuNum[p1]+=0.5*temp;
		      }
		    else
		      {
			OccuNum[p1]+=temp;
		      }
		  }
	    }
	}
    }

  for(int orb=0;orb<totalOrbitals;orb++)
    {
      if(orb>FermiSurface)
	{
	  OccuNum[orb]/=norm;
	}
      else
	{
	  OccuNum[orb]=(HFOrbs[orb].jj+1.)+OccuNum[orb]/norm;
	}
    }
}













//----------------------------------------------------------------
void MBPT::printHFOrbitals(int num) const
{
  if(num==-1)
    num=totalOrbitals;
  cout<<"num"<<setw(4)<<"\t"<<setw(4)<<"l"<<"\t"<<setw(4)<<"jj"<<"\t"<<setw(4)<<"tz"<<"\t"<<setw(10)<<"e"<<"\t"<<setw(4)<<"model space"<<endl;
  for(int i=0;i<num;i++)
    {
      cout<<i<<"\t"<<setw(4)<<HFOrbs[i].l<<"\t"<<setw(4)<<HFOrbs[i].jj<<"\t"<<setw(4)<<HFOrbs[i].tz<<"\t"<<setw(10)<<HFOrbs[i].e<<"\t"<<setw(4)<<ModelSpace[i]<<endl;
    }
}

void MBPT::printPOrbitals() const
{
  cout<<"num"<<setw(4)<<"\t"<<setw(4)<<"l"<<"\t"<<setw(4)<<"jj"<<"\t"<<setw(4)<<"tz"<<"\t"<<setw(10)<<"e"<<endl;
  for(int i=0;i<POrbitals.size();i++)
    {
      int orb=POrbitals[i];
      cout<<i<<"\t"<<setw(4)<<HFOrbs[orb].l<<"\t"<<setw(4)<<HFOrbs[orb].jj<<"\t"<<setw(4)<<HFOrbs[orb].tz<<"\t"<<setw(10)<<HFOrbs[orb].e<<endl;
    }
}


void MBPT::clearModelSpace()
{
  for(int i=0;i<ModelSpace.size();i++)
    ModelSpace[i]=false;
  POrbitals.clear();
}
void MBPT::setupPGroups()
{
  PGroups.clear();
  int NumOrbitals=POrbitals.size();
  GroupIndice.resize(NumOrbitals);
  OrbIndice.resize(NumOrbitals);
  for(int i=0;i<NumOrbitals;i++)
    {
      int j=0;
      for(j=0;j<PGroups.size();++j)
	{
	  if( PGroups[j].empty() ||
	      ( HFOrbs[POrbitals[i]].isInAGroupWith(HFOrbs[ PGroups[j][0] ]) )
	      )
	    {
	      PGroups[j].push_back(POrbitals[i]);
	      GroupIndice[i]=j;
	      OrbIndice[i]=PGroups[j].size()-1;
	      break;
	    }
	}
      if(j>=PGroups.size())
	{
	  PGroups.push_back(vector<int>() );
	  PGroups[PGroups.size()-1].push_back(POrbitals[i]);
	  GroupIndice[i]=PGroups.size()-1;
	  OrbIndice[i]=0;
	}
    } 
}

void MBPT::setupP2BStates()
{
  P2BStates.clear();
  for(int i=0;i<POrbitals.size();i++)
    for(int j=i;j<POrbitals.size();j++)
      {
	int a=POrbitals[i];
	int b=POrbitals[j];
	if(a>b) swap(a,b);
	int Jmin=abs(HFOrbs[a].jj - HFOrbs[b].jj)/2;
	int Jmax=(HFOrbs[a].jj + HFOrbs[b].jj)/2;
	for(int J=Jmin;J<=Jmax;J++)
	  {
	    int Tz=(HFOrbs[a].tz + HFOrbs[b].tz)/2;
	    int Par=(HFOrbs[a].l + HFOrbs[b].l)%2;
	    //if a==b, J must be even. |ab,J) = (-1)^{ja+jb-J+1} |ba,J)
	    if( (a==b) && J%2) continue;
	    P2BStates.push_back( TwoBodyStateType(a,b,Tz,Par,J) );
	  }
      }
  sort(P2BStates.begin(),P2BStates.end());
  //setup channels
  int NumTwoBodyStates=P2BStates.size();
  for(int i=0;i<NumTwoBodyStates;i++)
    {
      int LastChannel=P2BChannels.size()-1;
      //TwoBodyChannels is empty, so push_back a new channel.
      if(LastChannel<0)
	{
	  P2BChannels.push_back( vector<int>(1,i) );
	}
      else
	{
	  int indice=P2BChannels[LastChannel][0];
	  TwoBodyState & Temp=P2BStates[indice];
	  //the twobody state belong to LastChannel, so push_back the twobody state in this channel
	  if( P2BStates[i].isInAChannelWith(Temp) )
	    {
	      P2BChannels[LastChannel].push_back(i);
	    }
	  else//don't belong to, push_back a new channel
	    {
	      P2BChannels.push_back( vector<int>(1,i) );
	    }
	}
    }  
}

void MBPT::printP2BStates() const
{
  cout<<"num \t a \t b \t Tz \t Par \t J\n";
  for(int i=0;i<P2BStates.size();i++)
    {
      cout<<i<<"\t"<<P2BStates[i].a<<"\t"<<P2BStates[i].b<<"\t"<<P2BStates[i].Tz<<"\t"<<P2BStates[i].Par<<"\t"<<P2BStates[i].J<<endl;
    }
}
void MBPT::printP2BChannels() const
{
  cout<<"block\tJ\tTz\tPar\tdim\n";  
  for(int nb=0;nb<P2BChannels.size();nb++)
    {
      int J=P2BStates[P2BChannels[nb][0]].J;
      int Tz=P2BStates[P2BChannels[nb][0]].Tz;
      int Par=P2BStates[P2BChannels[nb][0]].Par;
      int dim=P2BChannels[nb].size();      
      cout<<nb<<"\t"<<J<<"\t"<<Tz<<"\t"<<Par<<"\t"<<dim<<endl;;
    }
}


MBPT::DataType MBPT::diagram0B(DataType w) 
{
  return HFSystem.HFenergy-CoreSumSPE+diagram0B_2p2h(w);//+diagram0B_4p2h(w)+diagram0B_2p4h(w)+diagram0B_3p3h(w);
}

MBPT::DataType MBPT::get0BQ(DataType w)
{
  return diagram0B(w);
}

void MBPT::get0BQDQ(DataType w,DataType & Q,DataType & DQ)
{
  Q=get0BQ(w);
  const double dw=1e-6;
  DQ=get0BQ(w+dw)-Q;
  DQ/=dw;
  //  DQ=diagram0B_2p2h_D1(w);
}

MBPT::DataType MBPT::get0BZ(DataType w)
{
  DataType DQ,Q;
  get0BQDQ(w,Q,DQ);
  return ( Q-DQ*w )/(1.-DQ);
}


void MBPT::solve0B()
{
  DataType w=0;
  DataType wold;
  double alpha=1.;
  do
    {
      wold=w;
      w=get0BZ(wold);
      w=alpha*w+(1-alpha)*wold;
      cout<<"0B:\t"<<w<<endl;
    }while(abs(w-wold)>1e-5);
  cout<<"0B:\t"<<w+CoreSumSPE<<endl;
  E0=w+CoreSumSPE;
  
  DataType Q,DQ;
  get0BQDQ(w,Q,DQ);
  rhoQ0B=-DQ;
  rhoQ0B=rhoQ0B/(1.+rhoQ0B);
}


MBPT::DataType MBPT::diagram1B(int a,int c,DataType w)
{
  return diagram1B_2p1h(a,c,w)+diagram1B_3p2h(a,c,w);
}
MBPT::DataType MBPT::get1BQ(int a,int c,DataType w)
{
  DataType temp=diagram1B(a,c,w);
  // if(a==c)
  //   temp+=diagram0B(w-HFOrbs[a].e);
  return temp;
}
void MBPT::get1BQ(int nb,Mat & Q,DataType w)
{
  int dim=PGroups[nb].size();
  Q.resize(dim,dim);
#pragma omp parallel for  
  for(int i=0;i<dim;i++)
    for(int j=i;j<dim;j++)
      {
	int a=PGroups[nb][i];
	int c=PGroups[nb][j];	
	Q(i,j)=get1BQ(a,c,w);
	Q(j,i)=Q(i,j);
      }
}
void MBPT::get1BQDQ(int nb,Mat & Q,Mat & DQ,DataType w)
{
  int dim=PGroups[nb].size();
  const double dw=1e-6;
  get1BQ(nb,Q,w);
  get1BQ(nb,DQ,w+dw);
  DQ-=Q;
  DQ/=dw;
}
void MBPT::get1BZ(int nb,Mat & Z,DataType w)
{
  int dim=PGroups[nb].size();
  Mat SP=Mat::Zero(dim,dim);
  for(int i=0;i<dim;i++)
    SP(i,i)=HFOrbs[PGroups[nb][i]].e;
  Mat I=Mat::Identity(dim,dim);
  Mat DQ,Q;
  get1BQDQ(nb,Q,DQ,w);
  Z=(I-DQ).inverse()*(Q-DQ*(w*I-SP));
}
void MBPT::solve1B()
{
  int blocks=PGroups.size();
  OneBodyVecs.resize(blocks);
  OneBodyVals.resize(blocks);
  P1BMat.resize(blocks);
  rhoQ1B.resize(blocks);
  for(int nb=0;nb<blocks;nb++)
    {
      int dim=PGroups[nb].size();
      Mat & h=P1BMat[nb];
      h=Mat::Zero(dim,dim);
      OneBodyVecs[nb]=Mat::Identity(dim,dim);
      OneBodyVals[nb].resize(dim);

      Mat SP=Mat::Zero(dim,dim);
      for(int i=0;i<dim;i++)
	{
	  SP(i,i)=HFOrbs[PGroups[nb][i]].e;
	  OneBodyVals[nb](i)=HFOrbs[PGroups[nb][i]].e;
	}
      Vec ValsOld;
      const double alpha=1.;
      int iter=0;
      do
	{
	  ValsOld=OneBodyVals[nb];
	  h=Mat::Zero(dim,dim);	  
	  for(int i=0;i<dim;i++)
	    {
	      Mat temp;
	      get1BZ(nb,temp,ValsOld(i));
	      temp+=SP;
	      Mat BiOrtho=OneBodyVecs[nb].inverse();
	      h+=temp*OneBodyVecs[nb].col(i)*BiOrtho.row(i);
	    }
	  ComplexEigenSolver<MatrixXcd> es(h);
	  // RealSchur<Mat> es(h);
	  // for(int i=0;i<dim;i++)
	  //   OneBodyVals[nb](i)=es.matrixT()(i,i);
	  // OneBodyVecs[nb]=es.matrixV();
	  OneBodyVals[nb]=alpha*es.eigenvalues()+(1-alpha)*ValsOld;
	  OneBodyVecs[nb]=es.eigenvectors();
	  quickSort(OneBodyVals[nb],OneBodyVecs[nb],0,dim-1);	  
	  cout<<nb<<"\t"<<"1B:\t"<<OneBodyVals[nb].transpose()<<endl;

	  // for(int i=0;i<dim;i++)
	  //   {
	  //     Mat mat1,mat2;
	  //     get1BZ(nb,mat1,OneBodyVals[nb](i));
	  //     get1BZ(nb,mat2,OneBodyVals[nb](i)+0.00001);
	  //     mat2-=mat1;
	  //     mat2/=0.00001;
	  //     cout<<(mat2*OneBodyVecs[nb].col(i)).transpose()<<endl;
	  //     cout<<endl;
	  //   }
	  // cout<<endl<<endl;

	}while( (ValsOld-OneBodyVals[nb]).norm()/dim>1e-6 );
      // cout<<"1B:\t"<<HFOrbs[PGroups[nb][0]].l<<"\t"<<HFOrbs[PGroups[nb][0]].jj<<"\t"<<HFOrbs[PGroups[nb][0]].tz<<"\t"<<( OneBodyVals[nb]+(CoreSumSPE-E0)*Vec::Ones(dim) ).transpose()<<endl;
      rhoQ1B[nb].resize(dim);
      Mat S,DS;
      for(int i=0;i<dim;i++)
	{
	  get1BQDQ(nb,S,DS,OneBodyVals[nb](i));
	  DataType rho=-OneBodyVecs[nb].col(i).transpose()*DS*OneBodyVecs[nb].col(i);
	  rhoQ1B[nb](i)=rho/(1.+rho);
	}
    }
}

void MBPT::printPGroups() const
{
  cout<<"block \t l \t jj \t tz \t dim\n";  
  for(int nb=0;nb<PGroups.size();nb++)
    {
      int jj=HFOrbs[PGroups[nb][0]].jj;
      int tz=HFOrbs[PGroups[nb][0]].tz;
      int l=HFOrbs[PGroups[nb][0]].l;
      int dim=PGroups[nb].size();      
      cout<<nb<<"\t"<<l<<"\t"<<jj<<"\t"<<tz<<"\t"<<dim<<endl;;
    }
}

void MBPT::printP1BMat() const
{
  for(int nb=0;nb<P1BMat.size();nb++)
    {
      int dim =P1BMat[nb].size();
      for(int i=0;i<dim;i++)
	for(int j=i;j<dim;j++)
	  {
	    //	    cout<<IndexInP[PGroups[nb][i]]<<"\t"<<IndexInP[PGroups[nb][j]]<<"\t"<<P1BMat[nb](i,j)+CoreSumSPE-E0<<endl;
	    cout<<IndexInP[PGroups[nb][i]]<<"\t"<<IndexInP[PGroups[nb][j]]<<"\t"<<P1BMat[nb](i,j)<<endl;
	  }
    }
}


MBPT::DataType MBPT::diagram2B(int a,int b,int c,int d,int J,DataType w)
{
  return get2B(a,b,c,d,J)+diagram2B_3p1h(a,b,c,d,J,w)+diagram2B_2p(a,b,c,d,J,w)+diagram2B_4p2h(a,b,c,d,J,w);
}
MBPT::DataType MBPT::diagram1B_to_2B(int a,int b,int c,int d,int J,DataType w)
{
  DataType temp=0;
  double phasecd= phase(J-(HFOrbs[c].jj+HFOrbs[d].jj)/2 + 1);  
  if( a==c && HFOrbs[b].isInAGroupWith(HFOrbs[d]) )
    {
      temp += diagram1B(b,d,w-HFOrbs[a].e);
    }
  if( b==d && HFOrbs[a].isInAGroupWith(HFOrbs[c]) )
    {
      temp += diagram1B(a,c,w-HFOrbs[b].e);
    }
  if( a==d && HFOrbs[b].isInAGroupWith(HFOrbs[c]) )
    {
      temp += phasecd * diagram1B(b,c,w-HFOrbs[a].e);
    }
  if( b==c && HFOrbs[a].isInAGroupWith(HFOrbs[d]) )
    {
      temp += phasecd * diagram1B(a,d,w-HFOrbs[b].e);
    }
  return temp;
}

MBPT::DataType MBPT::get2BQ(int bra,int ket,DataType w)
{
  int a=P2BStates[bra].a;
  int b=P2BStates[bra].b;
  int c=P2BStates[ket].a;
  int d=P2BStates[ket].b;
  int J=P2BStates[bra].J;//we don't check bra and ket has same quantum number here, they should be the same.

  double norm=1.;
  if(a==b) norm*=sqrt(2.);
  if(c==d) norm*=sqrt(2.);

  DataType temp=0;

  temp+=diagram2B(a,b,c,d,J,w);
  
  //1B diagrams to 2B
  temp+=diagram1B_to_2B(a,b,c,d,J,w);

  temp/=norm;
  // if(bra==ket)
  //   temp+=diagram0B(w-HFOrbs[d].e-HFOrbs[c].e);
  return temp;
}

void MBPT::get2BQ(int nb,Mat &Q, DataType w)
{
  int dim=P2BChannels[nb].size();
  Q.resize(dim,dim);
#pragma omp parallel for
  for(int i=0;i<dim;i++)
    for(int j=i;j<dim;j++)
      {
	int bra=P2BChannels[nb][i];
	int ket=P2BChannels[nb][j];
	Q(i,j)=get2BQ(bra,ket,w);
	Q(j,i)=Q(i,j);
      }
}

void MBPT::get2BQDQ(int nb,Mat & Q,Mat &DQ, DataType w)
{
  int dim=P2BChannels[nb].size();
  const double dw=1e-6;
  get2BQ(nb,Q,w);
  get2BQ(nb,DQ,w+dw);
  DQ-=Q;
  DQ/=dw;
}

void MBPT::get2BZ(int nb,Mat &Z, DataType w)
{
  int dim=P2BChannels[nb].size();  
  Mat SP=Mat::Zero(dim,dim);
  for(int i=0;i<dim;i++)
    {
      int ket=P2BChannels[nb][i];
      int a=P2BStates[ket].a;
      int b=P2BStates[ket].b;
      SP(i,i)=HFOrbs[a].e+HFOrbs[b].e;
    }
  Mat I=Mat::Identity(dim,dim);
  Mat DQ,Q;
  get2BQDQ(nb,Q,DQ,w);
  Z=(I-DQ).inverse()*(Q-DQ*(w*I-SP));  
}

void MBPT::solve2B()
{
  int blocks=P2BChannels.size();
  TwoBodyVecs.resize(blocks);
  TwoBodyVals.resize(blocks);
  P2BMat.resize(blocks);
  rhoQ2B.resize(blocks);

  for(int nb=0;nb<blocks;nb++)
    {
      int dim=P2BChannels[nb].size();
      Mat &h=P2BMat[nb];
      h=Mat::Zero(dim,dim);
      TwoBodyVecs[nb]=Mat::Identity(dim,dim);
      TwoBodyVals[nb].resize(dim);

      Mat SP=Mat::Zero(dim,dim);
      for(int i=0;i<dim;i++)
	{
	  int ket=P2BChannels[nb][i];
	  int a=P2BStates[ket].a;
	  int b=P2BStates[ket].b;
	  SP(i,i)=HFOrbs[a].e+HFOrbs[b].e;
	  TwoBodyVals[nb][i]=HFOrbs[a].e+HFOrbs[b].e;
	  //	  cout<<TwoBodyVals[nb][i]<<endl;
	}
      //init 
      // TwoBodyVals[nb][0]=10.7;
      // TwoBodyVals[nb][1]=12;
      // TwoBodyVals[nb][2]=15.;
      // for(int i=0;i<dim;i++)
      // 	{
      // 	  Mat Z;
      // 	  get2BZ(nb,Z,TwoBodyVals[nb][i]);
      // 	  Z+=SP;
      // 	  ComplexEigenSolver<Mat> esZ(Z);
      // 	  Vec vals=esZ.eigenvalues();
      // 	  Mat vecs=esZ.eigenvectors();
      // 	  quickSort(vals,vecs,0,dim-1);
      // 	  TwoBodyVecs[nb].col(i)=vecs.col(0);
      // 	  // cout<<Z<<endl;
      // 	  // cout<<vals<<endl;
      // 	  // cout<<vecs<<endl;
      // 	}
      
      Vec ValsOld;
      const double alpha=1.;
      int iter=0;
      do
	{
	  ValsOld=TwoBodyVals[nb];
	  h=Mat::Zero(dim,dim);	  
	  for(int i=0;i<dim;i++)
	    {
	      Mat temp;
	      get2BZ(nb,temp,ValsOld(i));
	      temp+=SP;
	      Mat BiOrtho=TwoBodyVecs[nb].inverse();
	      h+=temp*TwoBodyVecs[nb].col(i)*BiOrtho.row(i);
	    }

	  ComplexEigenSolver<MatrixXcd> es(h);
	  // RealSchur<Mat> es(h);
	  // for(int i=0;i<dim;i++)
	  //   TwoBodyVals[nb](i)=alpha*es.matrixT()(i,i)+(1-alpha)*ValsOld(i);
	  // TwoBodyVecs[nb]=es.matrixV();
	  TwoBodyVals[nb]=alpha*es.eigenvalues()+(1-alpha)*ValsOld;
	  TwoBodyVecs[nb]=es.eigenvectors();
	  quickSort(TwoBodyVals[nb],TwoBodyVecs[nb],0,dim-1);
	  cout<<nb<<"\t"<<"2B:\t"<<TwoBodyVals[nb].transpose()<<endl;

	  // for(int i=0;i<dim;i++)
	  //   {
	  //     Mat mat1,mat2;
	  //     get2BZ(nb,mat1,TwoBodyVals[nb](i));
	  //     get2BZ(nb,mat2,TwoBodyVals[nb](i)+0.00001);
	  //     mat2-=mat1;
	  //     mat2/=0.00001;
	  //     cout<<(mat2*TwoBodyVecs[nb].col(i)).transpose()<<endl;
	  //     cout<<endl;
	  //   }
	  // cout<<endl<<endl;
	  
	}while( (ValsOld-TwoBodyVals[nb]).norm()/dim>1e-5 );

      // for(int i=0;i<dim;i++)
      // 	{
      // 	  Mat temp;
      // 	  get2BZ(nb,temp,ValsOld(i));
      // 	  temp+=SP;
      // 	  ComplexEigenSolver<Mat> es(temp);
      // 	  cout<<ValsOld(i)<<"\t"<<es.eigenvalues().transpose()<<endl;
      // 	  cout<<endl;
      // 	}
      // Mat mat1,mat2;
      // get2BZ(nb,mat1,TwoBodyVals[nb](0));
      // get2BZ(nb,mat2,TwoBodyVals[nb](0)+0.00001);
      // mat2-=mat1;
      // mat2/=0.00001;
      // cout<<"2B:\t"<<P2BStates[P2BChannels[nb][0]].Par<<"\t"<<P2BStates[P2BChannels[nb][0]].J<<"\t"<<( TwoBodyVals[nb]+(CoreSumSPE-E0)*Vec::Ones(dim) ).transpose()<<endl;//"\t"<<mat2(0,0)<<endl;
      rhoQ2B[nb].resize(dim);
      Mat Q,DQ;
      for(int i=0;i<dim;i++)
	{
	  get2BQDQ(nb,Q,DQ,TwoBodyVals[nb](i));
	  DataType rho=-TwoBodyVecs[nb].col(i).transpose()*DQ*TwoBodyVecs[nb].col(i);
	  rhoQ2B[nb](i)=rho/(1.+rho);
	}
    }
}

MBPT::DataType MBPT::mat1B_to_2B(int a,int b,int c,int d,int J)
{
  DataType temp=0;
  double phasecd= phase(J-(HFOrbs[c].jj+HFOrbs[d].jj)/2 + 1);  
  if( a==c && HFOrbs[b].isInAGroupWith(HFOrbs[d]) )
    {
      int bra=OrbIndice[IndexInP[b]];
      int ket=OrbIndice[IndexInP[d]];
      temp += P1BMat[ GroupIndice[IndexInP[b]] ](bra,ket);
    }
  if( b==d && HFOrbs[a].isInAGroupWith(HFOrbs[c]) )
    {
      int bra=OrbIndice[IndexInP[a]];
      int ket=OrbIndice[IndexInP[c]];
      temp += P1BMat[ GroupIndice[IndexInP[a]] ](bra,ket);
    }
  if( a==d && HFOrbs[b].isInAGroupWith(HFOrbs[c]) )
    {
      int bra=OrbIndice[IndexInP[b]];
      int ket=OrbIndice[IndexInP[c]];
      temp += phasecd * P1BMat[ GroupIndice[IndexInP[b]] ](bra,ket);
    }
  if( b==c && HFOrbs[a].isInAGroupWith(HFOrbs[d]) )
    {
      int bra=OrbIndice[IndexInP[a]];
      int ket=OrbIndice[IndexInP[d]];
      temp += phasecd * P1BMat[ GroupIndice[IndexInP[a]] ](bra,ket);
    }
  double norm=1.;
  if(a==b) norm*=sqrt(2.);
  if(c==d) norm*=sqrt(2.);
  return temp/norm;
}

void MBPT::mat2Bminus1B()
{
  for(int nb=0;nb<P2BChannels.size();nb++)
    {
      int dim=P2BChannels[nb].size();
      for(int i=0;i<dim;i++)
	for(int j=0;j<dim;j++)
	  {
	    int bra=P2BChannels[nb][i];
	    int ket=P2BChannels[nb][j];
	    int a=P2BStates[bra].a;
	    int b=P2BStates[bra].b;
	    int c=P2BStates[ket].a;
	    int d=P2BStates[ket].b;
	    int J=P2BStates[bra].J;
	    P2BMat[nb](i,j)-=mat1B_to_2B(a,b,c,d,J);
	  }
    }
}
void MBPT::symmetrize()
{
  for(int nb=0;nb<P2BChannels.size();nb++)
    {
      int dim=P2BChannels[nb].size();
      Mat temp=P2BMat[nb];
      ComplexEigenSolver<MatrixXcd> es(temp);
      Mat U=es.eigenvectors();
      Mat T;
      InvSqrt(U*U.transpose(),T);
      T*=U;
      P2BMat[nb]=T * es.eigenvalues().asDiagonal() * T.transpose();
    }  
}

void MBPT::printResult() const
{
  // cout<<"****0B****\n";
  // cout<<E0<<endl;
  // cout<<rhoQ0B<<endl;
  // cout<<"========\n\n";    
  
  // cout<<"****1B****\n";
  // for(int nb=0;nb<PGroups.size();nb++)
  //   {
  //     int l=HFOrbs[PGroups[nb][0]].l;
  //     int jj=HFOrbs[PGroups[nb][0]].jj;
  //     int tz=HFOrbs[PGroups[nb][0]].tz;
  //     int dim=PGroups[nb].size();
  //     cout<<"--------\n";      
  //     cout<<l<<"\t"<<jj<<"\t"<<tz<<endl;
  //     cout<<( OneBodyVals[nb]+CoreSumSPE*Vec::Ones(dim) ).transpose()<<endl;
  //     cout<<rhoQ1B[nb].transpose()<<endl;
  //     cout<<"--------\n";      
  //   }
  // cout<<"========\n\n";
  
  cout<<"****2B****\n";
  for(int nb=0;nb<P2BChannels.size();nb++)
    {
      int J=P2BStates[P2BChannels[nb][0]].J;
      int Tz=P2BStates[P2BChannels[nb][0]].Tz;
      int Par=P2BStates[P2BChannels[nb][0]].Par;
      int dim=P2BChannels[nb].size();
      cout<<"--------\n";
      cout<<J<<"\t"<<Tz<<"\t"<<Par<<endl;
      //      cout<<(TwoBodyVals[nb]+CoreSumSPE*Vec::Ones(dim)).transpose()<<endl;
      cout<<TwoBodyVals[nb].transpose()<<endl;      
      cout<<rhoQ2B[nb].transpose()<<endl;      
      cout<<"--------\n";
    }
  cout<<"========\n\n";
}

void MBPT::print0BGraphData(ostream & fout,double start,double step,int n)
{
  for(int i=0;i<n;i++)
    {
      double w=start+i*step;
      //      fout<<w<<"\t"<<get0BZ(w)<<"\t"<<get0BQ(w)<<endl;
      fout<<w<<"\t"<<get0BZ(w).real()<<"\t"<<get0BQ(w).real()<<endl;
    }
}
void MBPT::print1BGraphData(ostream & fout,int nb,double start,double step,int n)
{
  int dim=PGroups[nb].size();
  Mat SP=Mat::Zero(dim,dim);
  for(int i=0;i<dim;i++)
    {
      SP(i,i)=HFOrbs[PGroups[nb][i]].e;
    }
  for(int i=0;i<n;i++)
    {
      double w=start+i*step;
      Mat Z,Q;
      get1BZ(nb,Z,w);
      get1BQ(nb,Q,w);	
      Z+=SP;
      Q+=SP;
      ComplexEigenSolver<Mat> esZ(Z);
      ComplexEigenSolver<Mat> esQ(Q);	
      Vec Zvals=esZ.eigenvalues();
      Vec Qvals=esQ.eigenvalues();
      quickSort(Zvals,0,dim-1);
      quickSort(Qvals,0,dim-1);      
      fout<<w<<"\t";
      for(int i=0;i<dim;i++)
	fout<<Zvals[i].real()<<"\t";
      for(int i=0;i<dim;i++)
	fout<<Qvals[i].real()<<"\t";
      fout<<endl;

    }
}

void MBPT::print2BGraphData(ostream & fout,int nb,double start,double step,int n)
{
  int J=P2BStates[P2BChannels[nb][0]].J;
  int Tz=P2BStates[P2BChannels[nb][0]].Tz;
  int Par=P2BStates[P2BChannels[nb][0]].Par;
  int dim=P2BChannels[nb].size();
  fout<<"----------------\n";
  fout<<"J\tTz\tPar\tdim\n";
  fout<<J<<"\t"<<Tz<<"\t"<<Par<<"\t"<<dim<<endl;
  fout<<"----------------\n";
  Mat SP=Mat::Zero(dim,dim);
  for(int i=0;i<dim;i++)
    {
      int ket=P2BChannels[nb][i];
      int a=P2BStates[ket].a;
      int b=P2BStates[ket].b;
      SP(i,i)=HFOrbs[a].e+HFOrbs[b].e;
    }
  for(int i=0;i<n;i++)
    {
      double w=start+i*step;
      Mat Z,Q;
      get2BZ(nb,Z,w);
      get2BQ(nb,Q,w);	
      Z+=SP;
      Q+=SP;
      ComplexEigenSolver<Mat> esZ(Z);
      ComplexEigenSolver<Mat> esQ(Q);	
      Vec Zvals=esZ.eigenvalues();
      Vec Qvals=esQ.eigenvalues();
      quickSort(Zvals,0,dim-1);
      quickSort(Qvals,0,dim-1);      
      
      fout<<w<<"\t";
      for(int i=0;i<dim;i++)
	fout<<Zvals[i].real()<<"\t";
      for(int i=0;i<dim;i++)
	fout<<Qvals[i].real()<<"\t";
      fout<<endl;
    }
  fout<<endl<<endl;
}




void MBPT::printP2BMat() const
{
  for(int nb=0;nb<P2BChannels.size();nb++)
    {
      int dim=P2BChannels[nb].size();
      for(int i=0;i<dim;i++)
	for(int j=i;j<dim;j++)
	  {
	    int bra=P2BChannels[nb][i];
	    int ket=P2BChannels[nb][j];
	    int a=P2BStates[bra].a;
	    int b=P2BStates[bra].b;
	    int c=P2BStates[ket].a;
	    int d=P2BStates[ket].b;
	    int J=P2BStates[bra].J;
	    //	    cout<<"\t"<<IndexInP[a]<<"\t"<<IndexInP[b]<<"\t"<<IndexInP[c]<<"\t"<<IndexInP[d]<<"\t"<<J<<"\t"<<P2BMat[nb](i,j)+( (i==j)?CoreSumSPE-E0:0 )<<endl;
	    cout<<"\t"<<IndexInP[a]+1<<"\t"<<IndexInP[b]+1<<"\t"<<IndexInP[c]+1<<"\t"<<IndexInP[d]+1<<"\t"<<J<<"\t"<<P2BMat[nb](i,j).real()<<"\t"<<P2BMat[nb](i,j).imag()<<endl;	    
	  }
    }
}










//------------------------diagrams

void MBPT::setIndice_hh(const TwoBodyStateType & State_hh,int channel,int indice_hh)
{
  int J=State_hh.J;
  int ha=State_hh.a;
  int hb=State_hh.b;
  int jjha=HFOrbs[ha].jj;
  int jjhb=HFOrbs[hb].jj;
  int hn_ab=J-abs(jjha-jjhb)/2;
    
  int hab=mapab_e(ha,hb,totalOrbitals);
  ChannelIndices_hhpp[hab][hn_ab] = channel;
  StateIndices_hhpp[hab][hn_ab] = indice_hh;
}

void MBPT::setIndice_pp(const TwoBodyStateType & State_pp,int channel,int indice_pp)
{
  int J=State_pp.J;
  int pa=State_pp.a;
  int pb=State_pp.b;
  int jjpa=HFOrbs[pa].jj;
  int jjpb=HFOrbs[pb].jj;
  int pn_ab=J-abs(jjpa-jjpb)/2;
    
  int pab=mapab_e(pa,pb,totalOrbitals);
  ChannelIndices_hhpp[pab][pn_ab] = channel;
  StateIndices_hhpp[pab][pn_ab] = indice_pp;
}

void MBPT::setIndice_hpm(const TwoBodyStateType & State_hpm,int channel,int indice_hpm)
{
  int h=State_hpm.a;
  int p=State_hpm.b;
  int J=State_hpm.J;
  int jjh=HFOrbs[h].jj;
  int jjp=HFOrbs[p].jj;
  int n_ab=J-abs(jjh-jjp)/2;
    
  ChannelIndices_hpm[h][p-FermiSurface-1][n_ab] = channel;
  StateIndices_hpm[h][p-FermiSurface-1][n_ab] = indice_hpm;
}
void MBPT::setIndice_phm(const TwoBodyStateType & State_phm,int channel,int indice_phm)
{
  int p=State_phm.a;
  int h=State_phm.b;
  int J=State_phm.J;
  int jjh=HFOrbs[h].jj;
  int jjp=HFOrbs[p].jj;
  int n_ab=J-abs(jjp-jjh)/2;
    
  ChannelIndices_phm[p-FermiSurface-1][h][n_ab] = channel;
  StateIndices_phm[p-FermiSurface-1][h][n_ab] = indice_phm;
}

void MBPT::setup2BStates()
{
  TwoBodyStates_hh.clear();
  TwoBodyStates_pp.clear();

  for(int a=0;a<=FermiSurface;a++)
    for(int b=a;b<=FermiSurface;b++)
      {
	const HFOrbital& A=HFOrbs[a];
	const HFOrbital & B=HFOrbs[b];
	int tz=(A.tz + B.tz)/2;
	int par=(A.l + B.l)%2;
	int jmin=abs(A.jj - B.jj)/2;
	int jmax=(A.jj + B.jj)/2;
	for(int j=jmin;j<=jmax;j++)
	  {
	    if(a==b && (j%2) ) continue;
	    TwoBodyStates_hh.push_back(TwoBodyStateType(a,b,tz,par,j));
	  }
      }
  for(int a=FermiSurface+1;a<totalOrbitals;a++)
    for(int b=a;b<totalOrbitals;b++)
      {
	const HFOrbital & A=HFOrbs[a];
	const HFOrbital & B=HFOrbs[b];
	int tz=(A.tz + B.tz)/2;
	int par=(A.l + B.l)%2;
	int jmin=abs(A.jj - B.jj)/2;
	int jmax=(A.jj + B.jj)/2;
	for(int j=jmin;j<=jmax;j++)
	  {
	    if(a==b && (j%2) ) continue;
	    TwoBodyStates_pp.push_back(TwoBodyStateType(a,b,tz,par,j));
	  }
      }
  sort(TwoBodyStates_hh.begin(),TwoBodyStates_hh.end());
  sort(TwoBodyStates_pp.begin(),TwoBodyStates_pp.end());


  
  TwoBodyStates_phm.clear();
  TwoBodyStates_hpm.clear();
  for(int a = FermiSurface+1; a <totalOrbitals; a++)  
    for(int b = 0; b <= FermiSurface; b++)
      {
	const HFOrbital & A=HFOrbs[a];
	const HFOrbital & B=HFOrbs[b];
	int tz=(A.tz - B.tz)/2;
	int par=abs(A.l - B.l)%2;
	int jmin=abs(A.jj - B.jj)/2;
	int jmax=(A.jj + B.jj)/2;
	for(int j=jmin;j<=jmax;j++)
	  TwoBodyStates_phm.push_back(TwoBodyStateType(a,b,tz,par,j));
      }

  for(int a = 0; a <=FermiSurface; a++)
    for(int b = FermiSurface+1; b <totalOrbitals; b++)
      {
	const HFOrbital & A=HFOrbs[a];
	const HFOrbital & B=HFOrbs[b];
	int tz=(A.tz - B.tz)/2;
	int par=abs(A.l - B.l)%2;
	int jmin=abs(A.jj - B.jj)/2;
	int jmax=(A.jj + B.jj)/2;
	for(int j=jmin;j<=jmax;j++)
	  TwoBodyStates_hpm.push_back(TwoBodyStateType(a,b,tz,par,j));
      }
  sort(TwoBodyStates_phm.begin(),TwoBodyStates_phm.end());
  sort(TwoBodyStates_hpm.begin(),TwoBodyStates_hpm.end());
}

void MBPT::setupVMat()
{
  V_hhhh.clear();
  V_hhpp.clear();
  V_pppp.clear();
  V_phmhpm.clear();
  V_phmphm.clear();
  e_hhpp.clear();

  ChannelIndices_hhpp.resize( (totalOrbitals*(totalOrbitals+1))/2 );
  StateIndices_hhpp.resize( (totalOrbitals*(totalOrbitals+1))/2 );
  for(int a=0;a<totalOrbitals;a++)
    for(int b=a;b<totalOrbitals;b++)
      {
	int jja=HFOrbs[a].jj;
	int jjb=HFOrbs[b].jj;
	int count= (jja+jjb)/2-abs(jja-jjb)/2 + 1;
	ChannelIndices_hhpp[mapab_e(a,b,totalOrbitals)].resize(count);
	StateIndices_hhpp[mapab_e(a,b,totalOrbitals)].resize(count);	
	for(int j=0;j<count;j++)
	  {
	    ChannelIndices_hhpp[mapab_e(a,b,totalOrbitals)][j]=-1;
	  }
      }
  
  //hhpp
  for (int i = 0, j = 0; i < TwoBodyStates_hh.size() && j < TwoBodyStates_pp.size(); i++)
    {
      TwoBodyStateType * State_hh = &TwoBodyStates_hh[i];
      TwoBodyStateType * State_pp = &TwoBodyStates_pp[j];
      while(!State_pp->isInAChannelWith(*State_hh) )
	{
	  j++;
	  if(j>=TwoBodyStates_pp.size() )
	    break;
	  State_pp = &TwoBodyStates_pp[j];
	}
      if( j >= TwoBodyStates_pp.size() )
	break;
      V_hhpp.push_back( Channel(State_hh->Tz,State_hh->Par,State_hh->J) );
      int channelNumber = V_hhpp.size()-1;

      V_hhpp[channelNumber].bra.push_back(i);
      V_hhpp[channelNumber].ket.push_back(j);
      setIndice_hh(*State_hh,channelNumber,V_hhpp[channelNumber].bra.size()-1);
      setIndice_pp(*State_pp,channelNumber,V_hhpp[channelNumber].ket.size()-1);

      if( i < TwoBodyStates_hh.size()-1 )
	{
	  TwoBodyStateType * Ph = & TwoBodyStates_hh[i+1];
	  while( (*Ph).isInAChannelWith(*State_hh) )
	    {
	      i++;
	      V_hhpp[channelNumber].bra.push_back(i);
	      setIndice_hh(*Ph,channelNumber,V_hhpp[channelNumber].bra.size()-1);
	      if( i >= TwoBodyStates_hh.size()-1) break;
	      State_hh = &TwoBodyStates_hh[i];
	      Ph = &TwoBodyStates_hh[i+1];
	    }
	}
      if(j < TwoBodyStates_pp.size()-1 )
	{
	  TwoBodyStateType * Pp = &TwoBodyStates_pp[j+1];
	  while( (*Pp).isInAChannelWith(*State_pp) )
	    {
	      j++;
	      V_hhpp[channelNumber].ket.push_back(j);
	      setIndice_pp(*Pp,channelNumber,V_hhpp[channelNumber].ket.size()-1);

	      if(j >= TwoBodyStates_pp.size()-1) break;
	      State_pp = &TwoBodyStates_pp[j];
	      Pp = &TwoBodyStates_pp[j+1];
	    }
	}
      j++;
    }

  //hhhh
  for(int i = 0; i < V_hhpp.size(); i++)
    {
      int tz = V_hhpp[i].Tz;
      int par = V_hhpp[i].Par;
      int j = V_hhpp[i].J;
      V_hhhh.push_back(Channel(tz,par,j));
      int size = V_hhpp[i].bra.size();
      V_hhhh[V_hhhh.size()-1].bra.resize(size);
      V_hhhh[V_hhhh.size()-1].ket.resize(size);
      for(int j = 0; j < size; j++)
	{
	  V_hhhh[V_hhhh.size()-1].bra[j] = V_hhpp[i].bra[j];
	  V_hhhh[V_hhhh.size()-1].ket[j] = V_hhpp[i].bra[j];
	}
    }

  //pppp
  for(int i = 0; i < V_hhpp.size(); i++)
    {
      int tz = V_hhpp[i].Tz;
      int par = V_hhpp[i].Par;
      int j = V_hhpp[i].J;
      V_pppp.push_back(Channel(tz,par,j));
      int size = V_hhpp[i].ket.size();
      V_pppp[V_pppp.size()-1].bra.resize(size);
      V_pppp[V_pppp.size()-1].ket.resize(size);
      for(int j = 0; j < size; j++)
	{
	  V_pppp[V_pppp.size()-1].bra[j] = V_hhpp[i].ket[j];
	  V_pppp[V_pppp.size()-1].ket[j] = V_hhpp[i].ket[j];
	}
    }
  
  //interaction for hhpp
  e_hhpp = V_hhpp;
#pragma omp parallel for
  for(int i = 0; i < V_hhpp.size(); i++)
    {
      int j=V_hhpp[i].J;
      int bra_dim = V_hhpp[i].bra.size();
      int ket_dim = V_hhpp[i].ket.size();
      V_hhpp[i].Mat.resize(bra_dim,ket_dim);
      e_hhpp[i].Mat.resize(bra_dim,ket_dim);
      for(int Ibra=0; Ibra < bra_dim; Ibra++)
	{
	  for(int Iket = 0; Iket < ket_dim; Iket++)
	    {
	      TwoBodyStateType & State_bra = TwoBodyStates_hh[ V_hhpp[i].bra[Ibra] ];
	      TwoBodyStateType & State_ket = TwoBodyStates_pp[ V_hhpp[i].ket[Iket] ];
	      int bra_a = State_bra.a;
	      int bra_b = State_bra.b;
	      int ket_a = State_ket.a;
	      int ket_b = State_ket.b;
	      double factor=1.0;
	      if(bra_a==bra_b)
		factor*=sqrt(0.5);
	      if(ket_a==ket_b)
		factor*=sqrt(0.5);
	      DataType V = get2B(bra_a,bra_b,ket_a,ket_b,j)*factor;
	      DataType de = HFOrbs[bra_a].e + HFOrbs[bra_b].e - HFOrbs[ket_a].e - HFOrbs[ket_b].e;
	      V_hhpp[i].Mat(Ibra,Iket) = V;
	      e_hhpp[i].Mat(Ibra,Iket) = de;
	    }
	}
    }

  //interaction for hhhh
#pragma omp parallel for  
  for(int i = 0; i < V_hhhh.size(); i++)
    {
      int j=V_hhhh[i].J;
      int dim = V_hhhh[i].bra.size();
      V_hhhh[i].Mat.resize(dim,dim);
      for(int Ibra=0; Ibra < dim; Ibra++)
	{
	  for(int Iket = Ibra; Iket < dim; Iket++)
	    {
	      TwoBodyStateType & State_bra = TwoBodyStates_hh[ V_hhhh[i].bra[Ibra] ];
	      TwoBodyStateType & State_ket = TwoBodyStates_hh[ V_hhhh[i].ket[Iket] ] ;
	      int bra_a = State_bra.a;
	      int bra_b = State_bra.b;
	      int ket_a = State_ket.a;
	      int ket_b = State_ket.b;
	      double factor=1.0;
	      if(bra_a==bra_b)
		factor*=sqrt(0.5);
	      if(ket_a==ket_b)
		factor*=sqrt(0.5);
	      V_hhhh[i].Mat(Ibra,Iket) = get2B(bra_a,bra_b,ket_a,ket_b,j)*factor;
	      V_hhhh[i].Mat(Iket,Ibra) = V_hhhh[i].Mat(Ibra,Iket);
	    }
	}
    }

  //interaction for pppp
#pragma omp parallel for  
  for(int i = 0; i < V_pppp.size(); i++)
    {
      int j=V_pppp[i].J;
      int dim = V_pppp[i].bra.size();
      V_pppp[i].Mat.resize(dim,dim);
      for(int Ibra = 0; Ibra < dim; Ibra++)
	{
	  for(int Iket = Ibra; Iket < dim; Iket++)
	    {
	      TwoBodyStateType & State_bra = TwoBodyStates_pp[ V_pppp[i].bra[Ibra] ];
	      TwoBodyStateType & State_ket = TwoBodyStates_pp[ V_pppp[i].ket[Iket] ];
	      int bra_a = State_bra.a;
	      int bra_b = State_bra.b;
	      int ket_a = State_ket.a;
	      int ket_b = State_ket.b;
	      double factor=1.0;
	      if(bra_a==bra_b)
		factor*=sqrt(0.5);
	      if(ket_a==ket_b)
		factor*=sqrt(0.5);
	      V_pppp[i].Mat(Ibra,Iket) = get2B(bra_a,bra_b,ket_a,ket_b,j)*factor;
	      V_pppp[i].Mat(Iket,Ibra) = V_pppp[i].Mat(Ibra,Iket);
	    }
	}
    }

  int numh=FermiSurface+1;
  int nump=totalOrbitals-FermiSurface;
  ChannelIndices_hpm.resize(numh);
  StateIndices_hpm.resize(numh);
  ChannelIndices_phm.resize(nump);
  StateIndices_phm.resize(nump);
  for(int i=0;i<=FermiSurface;i++)
    {
      ChannelIndices_hpm[i].resize(nump);
      StateIndices_hpm[i].resize(nump);
    }
  for(int j=FermiSurface+1;j<totalOrbitals;j++)
    {
      ChannelIndices_phm[j-FermiSurface-1].resize(numh);
      StateIndices_phm[j-FermiSurface-1].resize(numh);
    }
  for(int i=0;i<=FermiSurface;i++)
    for(int j=FermiSurface+1;j<totalOrbitals;j++)
      {
	int jja=HFOrbs[i].jj;
	int jjb=HFOrbs[j].jj;
	int count= (jja+jjb)/2-abs(jja-jjb)/2 + 1;
	ChannelIndices_hpm[i][j-FermiSurface-1].resize(count);
	StateIndices_hpm[i][j-FermiSurface-1].resize(count);
	ChannelIndices_phm[j-FermiSurface-1][i].resize(count);
	StateIndices_phm[j-FermiSurface-1][i].resize(count);	  
	for(int k=0;k<count;k++)
	  {
	    ChannelIndices_hpm[i][j-FermiSurface-1][k]=-1;
	    ChannelIndices_phm[j-FermiSurface-1][i][k]=-1;	      
	  }
      }
     
  
  //phmhpm
  for (int i = 0, j = 0; i < TwoBodyStates_phm.size() && j < TwoBodyStates_hpm.size(); i++)
    {
      TwoBodyStateType * State_phm = &TwoBodyStates_phm[i];      
      TwoBodyStateType * State_hpm = &TwoBodyStates_hpm[j];
      while( !State_hpm->isInAChannelWith(*State_phm) )
  	{
  	  j++;
  	  if(j >= TwoBodyStates_hpm.size() )
  	    break;
  	  State_hpm = &TwoBodyStates_hpm[j];
  	}
      if( j >= TwoBodyStates_hpm.size() )
  	break;

      V_phmhpm.push_back( Channel(State_phm->Tz,State_phm->Par,State_phm->J) );
      int channelNumber = V_phmhpm.size()-1;

      V_phmhpm[channelNumber].bra.push_back(i);
      V_phmhpm[channelNumber].ket.push_back(j);
      setIndice_phm(*State_phm,channelNumber,V_phmhpm[channelNumber].bra.size()-1);
      setIndice_hpm(*State_hpm,channelNumber,V_phmhpm[channelNumber].ket.size()-1);
      if( i < TwoBodyStates_phm.size()-1 )
	{
	  TwoBodyStateType * Pphm = &TwoBodyStates_phm[i+1];
	  while( (*Pphm).isInAChannelWith(*State_phm) )
	    {
	      i++;
	      V_phmhpm[channelNumber].bra.push_back(i);
	      setIndice_phm(*Pphm,channelNumber,V_phmhpm[channelNumber].bra.size()-1);
	      if( i >= TwoBodyStates_phm.size()-1) break;
	      State_phm = &TwoBodyStates_phm[i];
	      Pphm = &TwoBodyStates_phm[i+1];
	    }
	}
      if(j < TwoBodyStates_hpm.size()-1 )
	{
	  TwoBodyStateType * Phpm = &TwoBodyStates_hpm[j+1];
	  while( (*Phpm).isInAChannelWith(*State_hpm) )
	    {
	      j++;
	      V_phmhpm[channelNumber].ket.push_back(j);
	      setIndice_hpm(*Phpm,channelNumber,V_phmhpm[channelNumber].ket.size()-1);
	      if(j >= TwoBodyStates_hpm.size()-1) break;
	      State_hpm = &TwoBodyStates_hpm[j];
	      Phpm = &TwoBodyStates_hpm[j+1];
	    }
	}
      j++;
    }

  //phmphm
  for(int i = 0; i < V_phmhpm.size(); i++)
    {
      int tz = V_phmhpm[i].Tz;
      int par = V_phmhpm[i].Par;
      int J = V_phmhpm[i].J;
      V_phmphm.push_back(Channel(tz,par,J));
      int size = V_phmhpm[i].bra.size();
      V_phmphm[V_phmphm.size()-1].bra.resize(size);
      V_phmphm[V_phmphm.size()-1].ket.resize(size);
      for(int j = 0; j < size; j++)
	{
	  V_phmphm[V_phmphm.size()-1].bra[j] = V_phmhpm[i].bra[j];
	  V_phmphm[V_phmphm.size()-1].ket[j] = V_phmhpm[i].bra[j];
	}
    }

  //interaction for phmhpm  <p1 h1^-1 ,J | h2 p2^-1,J> define as (-1)^{jh2-jp2+J} cross_coupled1(p1,p2,h1,h2,J)
  e_phmhpm = V_phmhpm;
#pragma omp parallel for
  for(int i = 0; i < V_phmhpm.size(); i++)
    {
      int J=V_phmhpm[i].J;
      int bra_dim = V_phmhpm[i].bra.size();
      int ket_dim = V_phmhpm[i].ket.size();
      V_phmhpm[i].Mat.resize(bra_dim,ket_dim);
      e_phmhpm[i].Mat.resize(bra_dim,ket_dim);
      for(int Ibra=0; Ibra < bra_dim; Ibra++)
	{
	  for(int Iket = 0; Iket < ket_dim; Iket++)
	    {
	      TwoBodyStateType & State_bra = TwoBodyStates_phm[ V_phmhpm[i].bra[Ibra] ];
	      TwoBodyStateType & State_ket = TwoBodyStates_hpm[ V_phmhpm[i].ket[Iket] ];
	      int p1 = State_bra.a;//p1
	      int h1 = State_bra.b;//h1
	      int h2 = State_ket.a;//h2
	      int p2 = State_ket.b;//p2
	      DataType de = HFOrbs[h1].e + HFOrbs[h2].e - HFOrbs[p1].e - HFOrbs[p2].e;
	      int jjh2=HFOrbs[h2].jj;
	      int jjp2=HFOrbs[p2].jj;	      
	      V_phmhpm[i].Mat(Ibra,Iket) = cross_coupled1(p1,p2,h1,h2,J)*double(phase((jjh2-jjp2)/2+J));
	      e_phmhpm[i].Mat(Ibra,Iket) = de;
	    }
	}
    }
  
  //interaction for phmphm, <p1 h1^-1|p2 h2^-1 > define as (-1)^{jh2-jp2+J} cross_coupled1(p1,h2,h1,p2,J}
#pragma omp parallel for
  for(int i = 0; i < V_phmphm.size(); i++)
    {
      int J=V_phmphm[i].J;
      int dim = V_phmphm[i].bra.size();
      V_phmphm[i].Mat.resize(dim,dim);
      for(int Ibra=0; Ibra < dim; Ibra++)
	{
	  for(int Iket = Ibra; Iket < dim; Iket++)
	    {
	      TwoBodyStateType & State_bra = TwoBodyStates_phm[ V_phmphm[i].bra[Ibra] ];
	      TwoBodyStateType & State_ket = TwoBodyStates_phm[ V_phmphm[i].ket[Iket] ];
	      int p1 = State_bra.a;//p1
	      int h1 = State_bra.b;//h1
	      int p2 = State_ket.a;//p2
	      int h2 = State_ket.b;//h2
	      int jjh2=HFOrbs[h2].jj;
	      int jjp2=HFOrbs[p2].jj;
	      V_phmphm[i].Mat(Ibra,Iket) = cross_coupled1(p1,h2,h1,p2,J)*double(phase( (jjh2-jjp2)/2 + J));
	      V_phmphm[i].Mat(Iket,Ibra) = V_phmphm[i].Mat(Ibra,Iket);
	    }
	}
    }  
}

MBPT::DataType MBPT::getphmhpm(int p1,int h1,int h2,int p2,int J) const
{
  int jjp1=HFOrbs[p1].jj;
  int jjh1=HFOrbs[h1].jj;  
  int Jmin1=abs(jjp1-jjh1)/2;
  int channel1=ChannelIndices_phm[p1-FermiSurface-1][h1][J-Jmin1];
  int bra=StateIndices_phm[p1-FermiSurface-1][h1][J-Jmin1];  

  int jjh2=HFOrbs[h2].jj;
  int jjp2=HFOrbs[p2].jj;
  int Jmin2=abs(jjh2-jjp2)/2;
  int channel2=ChannelIndices_hpm[h2][p2-FermiSurface-1][J-Jmin2];
  int ket=StateIndices_hpm[h2][p2-FermiSurface-1][J-Jmin2];
  
  if(channel1!=channel2)
    return 0;
  if(channel1==-1)
    return 0;
  return V_phmhpm[channel1].Mat(bra,ket);
  //return cross_coupled1(h1,h2,p1,p2,J)*double(phase((jjh2-jjp2)/2+J));
}
MBPT::DataType MBPT::getphmphm(int p1,int h1,int p2,int h2,int J) const
{
  int jjp1=HFOrbs[p1].jj;
  int jjh1=HFOrbs[h1].jj;
  int Jmin1=abs(jjh1-jjp1)/2;
  int channel1=ChannelIndices_phm[p1-FermiSurface-1][h1][J-Jmin1];
  int bra=StateIndices_phm[p1-FermiSurface-1][h1][J-Jmin1];  

  int jjp2=HFOrbs[p2].jj;
  int jjh2=HFOrbs[h2].jj;
  int Jmin2=abs(jjh2-jjp2)/2;
  int channel2=ChannelIndices_phm[p2-FermiSurface-1][h2][J-Jmin2];
  int ket=StateIndices_phm[p2-FermiSurface-1][h2][J-Jmin2];
  
  if(channel1!=channel2)
    return 0;
  if(channel1==-1)
    return 0;
  return V_phmphm[channel1].Mat(bra,ket);
  //return phase((jjh2-jjp2)/2+J)*cross_coupled1(p1,h2,h1,p2,J);
}

MBPT::DataType MBPT::diagram0B_2p2h(DataType w)
{
  // DataType temp1=0,temp2=0,temp3=0;
  // for(int h1=0;h1<=FermiSurface;h1++)
  //   {
  //     for(int h2=h1;h2<=FermiSurface;h2++)
  // 	{

  // 	  int jj_min=abs(HFOrbs[h1].jj-HFOrbs[h2].jj);
  // 	  int jj_max=HFOrbs[h1].jj+HFOrbs[h2].jj;
  // 	  int Tz=HFOrbs[h1].tz+HFOrbs[h2].tz;
  // 	  int par=((HFOrbs[h1].l+HFOrbs[h2].l)%2);
  // 	  DataType e=HFOrbs[h1].e+HFOrbs[h2].e;
  // 	  for(int jj=jj_min;jj<=jj_max;jj+=2)
  // 	    {
  // 	      for(int p1=FermiSurface+1;p1<totalOrbitals;p1++)
  // 		for(int p2=p1;p2<totalOrbitals;p2++)
  // 		  {

  // 		    if((HFOrbs[p1].tz+HFOrbs[p2].tz)!=Tz) continue;
  // 		    if((HFOrbs[p1].l+HFOrbs[p2].l-par)%2) continue;
  // 		    if(abs(HFOrbs[p1].jj-HFOrbs[p2].jj)>jj) continue;
  // 		    if((HFOrbs[p1].jj+HFOrbs[p2].jj)<jj) continue;
  // 		    DataType temp=(jj+1.)*get2B(p1,p2,h1,h2,jj/2)*get2B(p1,p2,h1,h2,jj/2)/(w+e-HFOrbs[p1].e-HFOrbs[p2].e);
  // 		    if(p1==p2&&h1==h2)
  // 		      {
  // 			temp3+=temp;
  // 		      }
  // 		    else if(p1==p2||h1==h2)
  // 		      {
  // 			temp2+=temp;
  // 		      }
  // 		    else
  // 		      {
  // 			temp1+=temp;
  // 		      }
  // 		  }
  // 	    }
  // 	}
  //   }
  // temp2*=0.5;
  // temp3*=0.25;
  // temp1+=temp2;
  // temp1+=temp3;
  // return temp1;

  
  DataType temp=0;
  for(int channel=0;channel<V_hhpp.size();channel++)
    {
      Mat we=w*Mat::Ones(e_hhpp[channel].Mat.rows(),e_hhpp[channel].Mat.cols()) + e_hhpp[channel].Mat;
      temp+=( V_hhpp[channel].Mat.cwiseQuotient(we) * V_hhpp[channel].Mat.transpose() ).trace() * (2*V_hhpp[channel].J+1.);
    }
  return temp;
}

MBPT::DataType MBPT::diagram0B_2p2h_D1(DataType w)
{
  // DataType temp1=0,temp2=0,temp3=0;
  // for(int h1=0;h1<=FermiSurface;h1++)
  //   {
  //     for(int h2=h1;h2<=FermiSurface;h2++)
  // 	{

  // 	  int jj_min=abs(HFOrbs[h1].jj-HFOrbs[h2].jj);
  // 	  int jj_max=HFOrbs[h1].jj+HFOrbs[h2].jj;
  // 	  int Tz=HFOrbs[h1].tz+HFOrbs[h2].tz;
  // 	  int par=((HFOrbs[h1].l+HFOrbs[h2].l)%2);
  // 	  DataType e=HFOrbs[h1].e+HFOrbs[h2].e;
  // 	  for(int jj=jj_min;jj<=jj_max;jj+=2)
  // 	    {
  // 	      for(int p1=FermiSurface+1;p1<totalOrbitals;p1++)
  // 		for(int p2=p1;p2<totalOrbitals;p2++)
  // 		  {

  // 		    if((HFOrbs[p1].tz+HFOrbs[p2].tz)!=Tz) continue;
  // 		    if((HFOrbs[p1].l+HFOrbs[p2].l-par)%2) continue;
  // 		    if(abs(HFOrbs[p1].jj-HFOrbs[p2].jj)>jj) continue;
  // 		    if((HFOrbs[p1].jj+HFOrbs[p2].jj)<jj) continue;
  // 		    DataType temp=(jj+1.)*get2B(p1,p2,h1,h2,jj/2)*get2B(p1,p2,h1,h2,jj/2)/pow(w+e-HFOrbs[p1].e-HFOrbs[p2].e,2);
  // 		    if(p1==p2&&h1==h2)
  // 		      {
  // 			temp3+=temp;
  // 		      }
  // 		    else if(p1==p2||h1==h2)
  // 		      {
  // 			temp2+=temp;
  // 		      }
  // 		    else
  // 		      {
  // 			temp1+=temp;
  // 		      }
  // 		  }
  // 	    }
  // 	}
  //   }
  // temp2*=0.5;
  // temp3*=0.25;
  // temp1+=temp2;
  // temp1+=temp3;
  // return -temp1;
  DataType temp=0;
  for(int channel=0;channel<V_hhpp.size();channel++)
    {
      Mat we=w*Mat::Ones(e_hhpp[channel].Mat.rows(),e_hhpp[channel].Mat.cols()) + e_hhpp[channel].Mat;
      Mat we2=we.cwiseProduct(we);
      temp+=( V_hhpp[channel].Mat.cwiseQuotient(we2) * V_hhpp[channel].Mat.transpose() ).trace() * (2*V_hhpp[channel].J+1.);
    }
  return -temp;
}

MBPT::DataType MBPT::diagram0B_4p2h(DataType w)
{
  // DataType temp1=0,temp2=0,temp3=0,temp4=0;
  // for(int h1=0;h1<=FermiSurface;h1++)
  //   for(int h2=h1;h2<=FermiSurface;h2++)
  //     {
  // 	int jj_min=abs(HFOrbs[h1].jj-HFOrbs[h2].jj);
  // 	int jj_max=HFOrbs[h1].jj+HFOrbs[h2].jj;
  // 	int Tz=HFOrbs[h1].tz+HFOrbs[h2].tz;
  // 	int par=((HFOrbs[h1].l+HFOrbs[h2].l)%2);
  // 	DataType eh=HFOrbs[h1].e+HFOrbs[h2].e;
  // 	for(int jj=jj_min;jj<=jj_max;jj+=2)
  // 	  {
  // 	    for(int p1=FermiSurface+1;p1<totalOrbitals;p1++)
  // 	      {
  // 		for(int p2=p1;p2<totalOrbitals;p2++)
  // 		  {
  // 		    if((HFOrbs[p1].tz+HFOrbs[p2].tz)!=Tz) continue;
  // 		    if((HFOrbs[p1].l+HFOrbs[p2].l-par)%2) continue;
  // 		    if(abs(HFOrbs[p1].jj-HFOrbs[p2].jj)>jj) continue;
  // 		    if((HFOrbs[p1].jj+HFOrbs[p2].jj)<jj) continue;
  // 		    DataType ep12=HFOrbs[p1].e+HFOrbs[p2].e;
  // 		    for(int p3=FermiSurface+1;p3<totalOrbitals;p3++)
  // 		      {
  // 			for(int p4=p3;p4<totalOrbitals;p4++)
  // 			  {
  // 			    if((HFOrbs[p3].tz+HFOrbs[p4].tz)!=Tz) continue;
  // 			    if((HFOrbs[p3].l+HFOrbs[p4].l-par)%2) continue;
  // 			    if(abs(HFOrbs[p3].jj-HFOrbs[p4].jj)>jj) continue;
  // 			    if((HFOrbs[p3].jj+HFOrbs[p4].jj)<jj) continue;
  // 			    DataType ep34=HFOrbs[p3].e+HFOrbs[p4].e;
  // 			    DataType temp=(jj+1.)*get2B(h1,h2,p1,p2,jj/2)*get2B(p1,p2,p3,p4,jj/2)*get2B(p3,p4,h1,h2,jj/2)/(w+eh-ep12)/(w+eh-ep34);
  // 			    if(p1==p2&&p3==p4&&h1==h2)
  // 			      temp1+=temp;
  // 			    else if((p1==p2&&p3==p4)||(h1==h2&&p1==p2)||(p3==p4&&h1==h2))
  // 			      temp2+=temp;
  // 			    else if((p1==p2)||(h1==h2)||(p3==p4))
  // 			      temp3+=temp;
  // 			    else
  // 			      temp4+=temp;
  // 			  }

  // 		      }

  // 		  }

  // 	      }
  // 	  }
  //     }
  // temp1*=0.125;
  // temp2*=0.25;
  // temp3*=0.5;
  // temp1+=temp2;
  // temp1+=temp3;
  // temp1+=temp4;
  // return temp1;
  DataType temp=0;
  for(int channel=0;channel<V_hhpp.size();channel++)
    {
      Mat we=w*Mat::Ones(e_hhpp[channel].Mat.rows(),e_hhpp[channel].Mat.cols()) + e_hhpp[channel].Mat;
      temp+=( V_hhpp[channel].Mat.cwiseQuotient(we) * V_pppp[channel].Mat * V_hhpp[channel].Mat.cwiseQuotient(we).transpose() ).trace() * (2*V_hhpp[channel].J+1.);
    }
  return temp;  
}

MBPT::DataType MBPT::diagram0B_2p4h(DataType w)
{
  // DataType temp1=0,temp2=0,temp3=0,temp4=0;
  // for(int p1=FermiSurface+1;p1<totalOrbitals;p1++)
  //   for(int p2=p1;p2<=totalOrbitals;p2++)
  //     {
  // 	int jj_min=abs(HFOrbs[p1].jj-HFOrbs[p2].jj);
  // 	int jj_max=HFOrbs[p1].jj+HFOrbs[p2].jj;
  // 	int Tz=HFOrbs[p1].tz+HFOrbs[p2].tz;
  // 	int par=((HFOrbs[p1].l+HFOrbs[p2].l)%2);
  // 	DataType ep=HFOrbs[p1].e+HFOrbs[p2].e;
  // 	for(int jj=jj_min;jj<=jj_max;jj+=2)
  // 	  {
  // 	    for(int h1=0;h1<=FermiSurface;h1++)
  // 	      {
  // 		for(int h2=h1;h2<=FermiSurface;h2++)
  // 		  {
  // 		    if((HFOrbs[h1].tz+HFOrbs[h2].tz)!=Tz) continue;
  // 		    if((HFOrbs[h1].l+HFOrbs[h2].l-par)%2) continue;
  // 		    if(abs(HFOrbs[h1].jj-HFOrbs[h2].jj)>jj) continue;
  // 		    if((HFOrbs[h1].jj+HFOrbs[h2].jj)<jj) continue;
  // 		    DataType eh12=HFOrbs[h1].e+HFOrbs[h2].e;
  // 		    for(int h3=0;h3<=FermiSurface;h3++)
  // 		      {
  // 			for(int h4=h3;h4<=FermiSurface;h4++)
  // 			  {
  // 			    if((HFOrbs[h3].tz+HFOrbs[h4].tz)!=Tz) continue;
  // 			    if((HFOrbs[h3].l+HFOrbs[h4].l-par)%2) continue;
  // 			    if(abs(HFOrbs[h3].jj-HFOrbs[h4].jj)>jj) continue;
  // 			    if((HFOrbs[h3].jj+HFOrbs[h4].jj)<jj) continue;
  // 			    DataType eh34=HFOrbs[h3].e+HFOrbs[h4].e;
  // 			    DataType temp=(jj+1.)*get2B(h1,h2,p1,p2,jj/2)*get2B(p1,p2,h3,h4,jj/2)*get2B(h3,h4,h1,h2,jj/2)/(w+eh12-ep)/(w+eh34-ep);
  // 			    if(p1==p2&&h1==h2&&h3==h4)
  // 			      temp1+=temp;
  // 			    else if((p1==p2&&h1==h2)||(h3==h4&&p1==p2)||(h1==h2&&h3==h4))
  // 			      temp2+=temp;
  // 			    else if((p1==p2)||(h1==h2)||(h3==h4))
  // 			      temp3+=temp;
  // 			    else
  // 			      temp4+=temp;
  // 			  }

  // 		      }

  // 		  }

  // 	      }
  // 	  }
  //     }
  // temp1*=0.125;
  // temp2*=0.25;
  // temp3*=0.5;
  // temp1+=temp2;
  // temp1+=temp3;
  // temp1+=temp4;
  // return temp1;
  DataType temp=0;
  for(int channel=0;channel<V_hhpp.size();channel++)
    {
      Mat we=w*Mat::Ones(e_hhpp[channel].Mat.rows(),e_hhpp[channel].Mat.cols()) + e_hhpp[channel].Mat;
      temp+=( V_hhpp[channel].Mat.cwiseQuotient(we).transpose() * V_hhhh[channel].Mat * V_hhpp[channel].Mat.cwiseQuotient(we) ).trace() * (2*V_hhpp[channel].J+1.);
    }
  return temp;
}


MBPT::DataType MBPT::diagram0B_3p3h(DataType w)
{
  // DataType temp=0;
  // for(int h1=0;h1<=FermiSurface;h1++)
  //   for(int p1=FermiSurface+1;p1<totalOrbitals;p1++)
  //     {
  // 	int jj_min=abs(HFOrbs[p1].jj-HFOrbs[h1].jj);
  // 	int jj_max=HFOrbs[p1].jj+HFOrbs[h1].jj;
  // 	int Tz=HFOrbs[p1].tz-HFOrbs[h1].tz;
  // 	int par=((HFOrbs[p1].l+HFOrbs[h1].l)%2);
  // 	DataType ehp1=HFOrbs[h1].e-HFOrbs[p1].e;
  // 	for(int jj=jj_min;jj<=jj_max;jj+=2)
  // 	  {
  // 	    for(int h2=0;h2<=FermiSurface;h2++)
  // 	      for(int p2=FermiSurface+1;p2<totalOrbitals;p2++)
  // 		{
  // 		  if((HFOrbs[h2].tz-HFOrbs[p2].tz)!=Tz) continue;
  // 		  if((HFOrbs[p2].l+HFOrbs[h2].l-par)%2) continue;
  // 		  if(abs(HFOrbs[p2].jj-HFOrbs[h2].jj)>jj) continue;
  // 		  if((HFOrbs[p2].jj+HFOrbs[h2].jj)<jj) continue;
  // 		  for(int h3=h2;h3<=FermiSurface;h3++)
  // 		    for(int p3=(h2==h3)?p2:FermiSurface+1;p3<totalOrbitals;p3++)
  // 		      // for(int h3=0;h3<=FermiSurface;h3++)		      
  // 		      //   for(int p3=FermiSurface+1;p3<totalOrbitals;p3++)
  // 		      {
  // 			if((HFOrbs[h3].tz-HFOrbs[p3].tz)!=Tz) continue;
  // 			if((HFOrbs[p3].l+HFOrbs[h3].l-par)%2) continue;
  // 			if(abs(HFOrbs[p3].jj-HFOrbs[h3].jj)>jj) continue;
  // 			if((HFOrbs[p3].jj+HFOrbs[h3].jj)<jj) continue;
  // 			int count;
  // 			count=2;
  // 			if(h2==h3&&p2==p3)
  // 			  count=1;
  // 			temp+=count*phase((HFOrbs[h2].jj-HFOrbs[p2].jj+jj)/2)*(jj+1.)*cross_coupled1(h1,h3,p1,p3,jj/2)*cross_coupled1(p2,h3,h2,p3,jj/2)*cross_coupled1(h1,h2,p1,p2,jj/2)/((w+HFOrbs[h1].e+HFOrbs[h2].e-HFOrbs[p1].e-HFOrbs[p2].e)*(w+HFOrbs[h1].e+HFOrbs[h3].e-HFOrbs[p1].e-HFOrbs[p3].e));
  // 		      }
  // 		}
  // 	  }
  //     }
  // return temp;
  DataType temp=0;
  for(int channel=0;channel<V_phmhpm.size();channel++)
    {
      Mat we=w*Mat::Ones(e_phmhpm[channel].Mat.rows(),e_phmhpm[channel].Mat.cols()) + e_phmhpm[channel].Mat;
      temp+=( V_phmhpm[channel].Mat.cwiseQuotient(we) * V_phmhpm[channel].Mat.cwiseQuotient(we).transpose() * V_phmphm[channel].Mat).trace() * (2*V_phmhpm[channel].J+1.);
    }
  return temp;  
}






MBPT::DataType MBPT::diagram1B_2p1h(int a,int c,DataType w)
{
  DataType temp=0;
  for(int h=0;h<=FermiSurface;h++)
    {
      int Jmin=abs(HFOrbs[h].jj-HFOrbs[a].jj)/2;
      int Jmax=(HFOrbs[h].jj+HFOrbs[a].jj)/2;
      int Tz=HFOrbs[h].tz+HFOrbs[a].tz;
      int par=((HFOrbs[h].l+HFOrbs[a].l)%2);
      for(int J=Jmin;J<=Jmax;J++)
	{
	  double factr=(J*2+1.)/(HFOrbs[a].jj+1.);
	  for(int p1=FermiSurface+1;p1<totalOrbitals;p1++)
	    for(int p2=p1;p2<totalOrbitals;p2++)
	      {
		if((HFOrbs[p1].tz+HFOrbs[p2].tz)!=Tz) continue;
		if((HFOrbs[p1].l+HFOrbs[p2].l-par)%2) continue;
		if(abs(HFOrbs[p1].jj-HFOrbs[p2].jj)>J*2) continue;
		if((HFOrbs[p1].jj+HFOrbs[p2].jj)<J*2) continue;
		double fac=1.0;
		if(p2==p1) fac=0.5;
		DataType de=w+HFOrbs[h].e-HFOrbs[p1].e-HFOrbs[p2].e;
		temp+=factr*fac*get2B(p1,p2,c,h,J)*get2B(a,h,p1,p2,J)/de;
	      }
	}
    }
  return temp;
}

MBPT::DataType MBPT::diagram1B_3p2h(int a,int c,DataType w)
{
  DataType temp=0;
  for(int p=FermiSurface+1;p<totalOrbitals;p++)
    {
      int Jmin=abs(HFOrbs[p].jj-HFOrbs[a].jj)/2;
      int Jmax=(HFOrbs[p].jj+HFOrbs[a].jj)/2;
      int Tz=HFOrbs[p].tz+HFOrbs[a].tz;
      int par=((HFOrbs[p].l+HFOrbs[a].l)%2);
      for(int J=Jmin;J<=Jmax;J++)
	{
	  double factr=(J*2+1.)/(HFOrbs[a].jj+1.);
	  for(int h1=0;h1<=FermiSurface;h1++)
	    for(int h2=h1;h2<=FermiSurface;h2++)
	      {
		if((HFOrbs[h1].tz+HFOrbs[h2].tz)!=Tz) continue;
		if((HFOrbs[h1].l+HFOrbs[h2].l-par)%2) continue;
		if(abs(HFOrbs[h1].jj-HFOrbs[h2].jj)>J*2) continue;
		if((HFOrbs[h1].jj+HFOrbs[h2].jj)<J*2) continue;
		double fac=1.0;
		if(h2==h1) fac=0.5;
		DataType de=w-HFOrbs[p].e-HFOrbs[a].e-HFOrbs[c].e+HFOrbs[h1].e+HFOrbs[h2].e;
		temp+=factr*fac*get2B(h1,h2,c,p,J)*get2B(a,p,h1,h2,J)/de;
	      }
	}
    }
  return temp;
}


MBPT::DataType MBPT::core_polarization(int a,int b,int c,int d,int J,DataType w)
{
  DataType temp=0;
  int jja=HFOrbs[a].jj;
  int jjb=HFOrbs[b].jj;
  int jjc=HFOrbs[c].jj;
  int jjd=HFOrbs[d].jj;  
  int phase1=phase((HFOrbs[b].jj+HFOrbs[c].jj)/2 +J);
  int Jminac=abs(HFOrbs[a].jj-HFOrbs[c].jj)/2;
  int Jmaxac=(HFOrbs[a].jj+HFOrbs[c].jj)/2;
  int Jminbd=abs(HFOrbs[b].jj-HFOrbs[d].jj)/2;
  int Jmaxbd=(HFOrbs[b].jj+HFOrbs[d].jj)/2;
  int Jmin=max(Jminac,Jminbd);
  int Jmax=min(Jmaxac,Jmaxbd);
  int Tz=HFOrbs[a].tz-HFOrbs[c].tz;
  int par=HFOrbs[a].l+HFOrbs[c].l;
  for(int Jph=Jmin;Jph<=Jmax;Jph++)
    {
      double sixj=gsl_sf_coupling_6j(jja,jjc,2*Jph,jjd,jjb,2*J);
      for(int h=0;h<=FermiSurface;h++)
	for(int p=FermiSurface+1;p<totalOrbitals;p++)
	  {
	    int jjh=HFOrbs[h].jj;
	    int jjp=HFOrbs[p].jj;
	    if((HFOrbs[p].tz-HFOrbs[h].tz)!=Tz) continue;
	    if((HFOrbs[p].l+HFOrbs[h].l-par)%2) continue;
	    if(abs(jjp-jjh)>Jph*2) continue;
	    if((jjp+jjh)<Jph*2) continue;
	    double phase2=phase((jjh-jjp)/2+Jph+Jph);
	    DataType de=w+HFOrbs[h].e-HFOrbs[p].e-HFOrbs[c].e-HFOrbs[b].e;
	    temp+=cross_coupled1(a,h,c,p,Jph)*cross_coupled1(p,b,h,d,Jph)*(2*Jph+1.)*(sixj*phase1*phase2)/de;
	  }
    }
  return temp;
}

MBPT::DataType MBPT::diagram2B_3p1h(int a,int b,int c,int d,int J,DataType w)
{
  double phaseab=phase( (HFOrbs[a].jj+HFOrbs[b].jj)/2-J+1 );
  double phasecd=phase( (HFOrbs[c].jj+HFOrbs[d].jj)/2-J+1 );
  DataType temp1,temp2,temp3,temp4;
  temp1=core_polarization(a,b,c,d,J,w);
  if(c==d)
    temp2=temp1*phasecd;
  else
    temp2=core_polarization(a,b,d,c,J,w)*phasecd;
  if(a==b)
    {
      temp3=temp1*phaseab;
      temp4=temp2*phaseab;
    }
  else
    {
      temp3=core_polarization(b,a,c,d,J,w)*phaseab;
      temp4=core_polarization(b,a,d,c,J,w)*phaseab*phasecd;
    }
  return temp1+temp2+temp3+temp4;
}

MBPT::DataType MBPT::diagram2B_2p(int a,int b,int c,int d,int J,DataType w)
{
  DataType temp=0;
  int Tz=HFOrbs[a].tz+HFOrbs[b].tz;
  int par=HFOrbs[a].l+HFOrbs[b].l;  
  for(int p1=FermiSurface+1;p1<totalOrbitals;p1++)
    for(int p2=p1;p2<totalOrbitals;p2++)
      {
	if( ModelSpace[p1]&&ModelSpace[p2] ) continue;
	if( (HFOrbs[p1].tz+HFOrbs[p2].tz)!=Tz ) continue;
	if((HFOrbs[p1].l+HFOrbs[p2].l-par)%2) continue;
	if( (HFOrbs[p1].jj+HFOrbs[p2].jj) < 2*J) continue;
	if( abs(HFOrbs[p1].jj-HFOrbs[p2].jj) > 2*J) continue;	
	DataType de=w-HFOrbs[p1].e-HFOrbs[p2].e;
	double fac=1.;
	if(p1==p2) fac=0.5;
	temp+=fac*get2B(a,b,p1,p2,J)*get2B(c,d,p1,p2,J)/de;
      }
  return temp;
}

MBPT::DataType MBPT::diagram2B_4p2h(int a,int b,int c,int d,int J,DataType w)
{
  DataType temp=0;
  int Tz=HFOrbs[a].tz+HFOrbs[b].tz;
  int par=HFOrbs[a].l+HFOrbs[b].l;  
  for(int h1=0;h1<=FermiSurface;h1++)
    for(int h2=h1;h2<=FermiSurface;h2++)
      {
	if( (HFOrbs[h1].tz+HFOrbs[h2].tz)!=Tz ) continue;
	if((HFOrbs[h1].l+HFOrbs[h2].l-par)%2) continue;
	if( (HFOrbs[h1].jj+HFOrbs[h2].jj) < 2*J) continue;
	if( abs(HFOrbs[h1].jj-HFOrbs[h2].jj) > 2*J) continue;	
	double fac=1.;
	if(h1==h2) fac=0.5;
	DataType de=w+HFOrbs[h1].e+HFOrbs[h2].e-HFOrbs[a].e-HFOrbs[b].e-HFOrbs[c].e-HFOrbs[d].e;
	temp+=fac*get2B(h1,h2,c,d,J)*get2B(h1,h2,a,b,J)/de;
      }
  return temp;
}
